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Abstract 

We use a semi-numerical method to find the position and period of periodic orbits in a bisymmetrical potential, made 
up of a two dimensional harmonic oscillator, with an additional term of a Plummer potential, in a number of resonant 
cases. The results are compared with the outcomes obtained by the numerical integration of the equations of motion and 
the agreement is good. This indicates that the semi-numerical method gives general and reliable results. Comparison 
with other methods of locating periodic orbits is also made. 

Keywords: Galaxies: kinematics and dynamics; periodic orbits 
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\1. Introduction 

Dynamical systems made up of harmonic oscillators 
have been extensively used over the last five decades in or- 
der to describe motion near an equilibrium point (see e.g 
Hcnon and Heiles, 1964; Giorgilli and Galgani, 1978; Saito 
and Ichimura, 1979; Caranicolas, 1993; 1994a; Caranico- 
las and Karanis, 1998; 1999; Arribas at al., 2006). In 
order to study these systems, scientists have used numeri- 
cal (Karanis and Vozikis, 2008) or analytical methods (see 
Caranicolas and Barbanis, 1982; Deprit, 1991; Deprit and 
Elipe, 1991). 

Astronomers frequently use potentials made up of har- 
monic oscillators, in order to study local motion in galax- 
ies. Of particular interest are the bisymmetrical potentials, 
as these systems have been used in order to describe mo- 
tion near the center of an elliptical galaxy. A large part of 
the above studies have been devoted to locate the position 
and find the period of the periodic orbits, as these orbits 
represent the backbone of the whole set of orbits. 

In an earlier paper Caranicolas and Innanen (1992) 
(hereafter Paper I), we studied the periodic motion in po- 
tentials of the form 



1 



Vi = - {uj\x 2 + uj\y 2 ) + h.o.t, 



(1) 



in the case where uj\ = W2, that is when the frequencies 
of oscillations along the x and y axis, respectively, are 
equal. This is called a perturbed elliptic oscillator. In 
paper I, we used a combination of numerical and analytical 
calculations in order to find the position and the period of 
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the periodic orbits for different h.o.t (higher order terms), 
that is for different perturbation functions. 
In this article we shall use the potential 



V(x,y) 



(cufx 2 



y/ x 2 + y 2 + , 



(2) 



where \x and a are parameters. Potential (2) consists of 
a two dimensional harmonic oscillator with an additional 
Plummer sphere and can be used to describe motion in the 
central parts of a galaxy with a central bulge and a dense 
nucleus. The equations of motion for a test particle with 
a unit mass are 



(x 2 + y 2 + a 2 ) 



3/2 



y 



(x 2 + y 2 + a 2 ) 



3/2 



y 



(3) 



where the dot indicates derivative with respect to the time. 
The corresponding Hamiltonian is written as 



H 



1 



/ 2 1 2 1 2 2 1 2 2\ 
(p x +p y +UJ 1 X +uj 2 y ) 



x 2 + y 2 + . 



h, 



(4) 



where p x and p y are the momenta per unit mass conju- 
gate to x and y, while h > is the numerical value of the 
Hamiltonian, which is conserved. Our aim is to find the 
position of periodic orbits and the corresponding periods 
in the resonant cases, using a semi-numerical procedure 
and also to compare the results with the outcomes given 
by the numerical integration of the equations of motion. 
We shall study the resonant cases Wi:w2 = 1:1) wi:w2 = 
2:1, u)\:u)2 = 2:3 and U3\:U2 = 4:3. In order to keep thing 
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simple, we use the values: /i = 0.001 and a = 0.25, while 
the value of the energy h, will be treated as a parameter. 
For the numerical integration of the equations of motion, a 
Bulirsh-Stocr numerical integration routine in double pre- 
cision is used. The accuracy of the calculations is checked 
by the constancy of the energy integral (4), which is con- 
served up to the twelfth significant figure. 

The paper is organized as follows: In Section 2 we find 
the position and the periods of the periodic orbits starting 
perpendicularly from the x axis. In Section 3 we study 
periodic orbits going through the origin. A discussion and 
the conclusions of this work are presented in Section 4. 

2. Periodic orbits starting perpendicularly from the 

x axis 

We shall start for case when uj\ = u 2 = u, that is 
for the 1:1 resonant case. In this case the potential (2) is 
axially symmetric and in polar coordinates (r, <p) takes the 
form 



2„2 



V(r) = -ufr 



\Jr 2 + a 2 



(5) 



As the potential is axially symmetric the test particle's 
angular momentum L is an integral of motion. In this 
case the system has a circular periodic orbit of radius r c . 
Elementary calculations show that the angular momentum 
of this circular orbit is given by the equation 



dV 



dr 



while the corresponding value of the energy is 
V(r). 



1 T 2 

h - 

2 rj. 



(6) 



(7) 



The period of the circular periodic orbit can be found, 
if we observe that as the orbit is a circle with a radius 
r c the two frequencies uj\ a and uj 2a in equations (3) are 
always equal with a value 



Win — ^>2a — — 



(J, 



(r 2 + a 2 ) 



3/2 



1/2 



Thus the period of the circular orbit is given by 
2tt 



LU 2 + 



(r=+a 2 ) 



2 -,3/2 



1/2- 



(8) 



(9) 



We now proceed to find the position of the periodic 
orbits in the case when cj\ — 2uj 2 , that is in the resonance 
case 2:1. Let x s0 be the starting position of the periodic 
orbit from the x axis. Remember that yo = p x o = 0, while 
the value of p y0 is always found from the energy integral 



(4) . Our numerical calculations show that the above value 
of p y o is in a good agreement with the value 



2^0 
UJ 2 ' 



(10) 



If we set x s q = x, y = p x = 0,p y = x s q/lo 2 , in the harmonic 
part of the Hamiltonian (4), that is when fi = and solve 
the corresponding equation for x s , we find 



1 + bj{ljJ 2 



1/2 



(11) 



which gives a first value for the starting point of the pe- 
riodic orbit. Numerical experiments indicate that a good 
approximation for the frequency of the 2:1 periodic orbits 
is 



Wis = 



tol + ■ 



fj, 



2(x 



\3/2 



1/2 



(12) 



Setting this value of the frequency in equation (10) we find 

(13) 



Wis 



where x s is the new value of the starting point of the 2:1 
periodic orbit. We now set x s = x,y = p x = 7 p ys = 
x s /w 2si w 2s = w\s/ 2 in Hamiltonian (4) and get 



1 



2 {Wls 



w 



2.s 



^fx 2 s +a 2 



(14) 



If we solve equation (14) for x s we obtain the starting point 
of the 2:1 periodic orbit. The period of the 2:1 periodic 
orbit is given by 



4-7T 
Wl s 



(15) 



Table 1 gives the positions and periods of the 2:1 pe- 
riodic orbits for several values of the energy h. Subscript 
n indicates values found by numerical integration, while 
subscript s indicates values found using equations (14) and 
(15). One can see that the differences in the position of pe- 
riodic orbits vary from to about 9%, while in all cases the 
differences in the period are less than 1%. The values of the 
parameters are: uj\ = 0.8, uj 2 = 0A,fi = 0.001, a = 0.25. 
Figure 1 shows a 2:1 periodic orbit, for the above values 
of parameters when h = 0.05. The initial conditions are: 
xq = x n = 0.1295, yo = p x o = 0, while the value of p y0 is 
found from the energy integral (4). The outermost curve 
is the curve of zero velocity. 

Next we come to study the 2:3 resonance, that is when 
(Ji : uj 2 = 2:3. Here the numerical calculations suggest 
that a good approximation for the starting point of the 
periodic orbit is given by the formula 



1 



UJ 2 



h, 



(16) 
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Table 1: Positions and periods of the 2:1 resonant periodic orbits 
intersecting the x axis perpendicularly. The values of the parameters 
are lo\ = 0.8 and LU2 = 0.4. Subscript s indicates values derived by 
semi-numerical methods, while subscript n shows results obtained 
by numerical integration. 



h 


X s 




T s 


T n 


0.049 


0.1254 


0.1254 


15.4270 


15.5148 


0.050 


0.1265 


0.1295 


15.4285 


15.5183 


0.051 


0.1276 


0.1335 


15.4300 


15.5217 


0.052 


0.1288 


0.1373 


15.4316 


15.5249 


0.053 


0.1300 


0.1410 


15.4331 


15.5281 


0.054 


0.1311 


0.1446 


15.4345 


15.5312 
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Figure 1: A 2:1 resonant periodic orbit intersecting the x axis per- 
pendicularly. The values of parameters are as in Table 1. The value 
of energy is h = 0.05. Initial conditions are: xq — Xn = 0.1295, yo = 
PxO = 0, while the value of p y o is found from the energy integral. 



while the value of the corresponding period is 

m 4-7T 

T s = , 

Wis 

where 

1/2 



Wl s = 



[1 



(x 2 s + a 2 ) 3/2 



(17) 



(18) 



Table 2 is similar to Table 1, but for the 2:3 periodic 
orbits. We see that the agreement between the results 
given by numerical integration and the outcomes from the 
semi-numerical formulas (16) and (17) is very good. The 
values of the parameters are: u>i — 0.4,^2 = 0.6, /i = 
0.001, a = 0.25. Figure 2 shows a 2:3 periodic orbit, for 
the above values of parameters when h = 0.08. Initial 



Table 2: Similar to Table 1 but for the 2:3 resonant periodic orbits. 
The values of the parameters arc: ui = 0.4, u>2 = 0.6. 



h 


x s 


x n 


T s 


T n 


0.065 


0.7806 


0.7832 


31.2391 


31.1167 


0.070 


0.8100 


0.8096 


31.2560 


31.1379 


0.075 


0.8385 


0.8334 


31.2704 


31.1578 


0.080 


0.8660 


0.8562 


31.2827 


31.1755 


0.085 


0.8926 


0.8789 


31.2934 


31.1911 


0.090 


0.9185 


0.9002 


31.3027 


31.2053 
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Figure 2: A 2:3 resonant periodic orbit intersecting the x axis per- 
pendicularly. The values of parameters are as in Table 2. The value 
of energy is h = 0.08. Initial conditions are: x$ = x n = 0.8562, = 
PxO = 0, while the value of p y o is found from the energy integral. 



conditions are: x$ = x n = 0.8562, j/o = PxO = 0, while the 
value of p y o is found from the energy integral (4). 

Finally we shall study the case when u>\ : u>2 = 4:3. In 
this case our numerical experiments show that the value 
of the starting point of the periodic orbit is well approxi- 
mated if we use the formula (16). The value of the corre- 
sponding period is now given by the formula 



8?r 

Wis 



where 



Wis = 



tot + 



2{xl + a 2 ) 



3/2 



1/2 



(19) 



(20) 



In Table 3 we make a comparison between the results 
given by the numerical integration and those using formu- 
las (16) and (19) for several values of the energy h. We see 
that the agreement is again very good. In all numerical 
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Tabic 3: Similar to Table 1 but for the 4:3 resonant periodic orbits. 
The values of the parameters arc: uii = 0.8,^2 = 0.6. 



h 


x s 






T n 


0.35 


0.6404 


0.6401 


31.3782 


31.3684 


0.40 


0.6846 


0.6847 


31.3843 


31.3762 


0.45 


0.7261 


0.7266 


31.3889 


31.3820 


0.50 


0.7654 


0.7656 


31.3924 


31.3865 


0.55 


0.8028 


0.8027 


31.3953 


31.3901 


0.60 


0.8385 


0.8389 


31.3976 


31.3931 



Table 4: Periods of the 1:1 periodic orbits going through the origin. 
The values of the parameters arc: uii = LU2 = 0.4. 



h 


T s 




0.02 


15.2813 


15.2839 


0.04 


15.4961 


15.4866 


0.06 


15.5770 


15.5644 


0.08 


15.6171 


15.6043 


0.10 


15.6403 


15.6280 


0.12 


15.6551 


15.6436 
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Figure 3: A 4:3 resonant periodic orbit intersecting the x axis per- 
pendicularly. The values of parameters arc as in Table 3. The value 
of energy is h = 0.40. Initial conditions are: xq — Xn = 0.6847, yo = 
PxO = 0, while the value of p y g is found from the energy integral. 



calculations we use the values: wi = 0.8, L02 = 0.6, /j = 
0.001, a = 0.25. Figure 3 shows a 4:3 periodic orbit, for 
the above values of parameters, when h = 0.40. Initial 
conditions arc: xq = x n = 0.6847, yo = p X Q = 0, while the 
value of p y o is found from the energy integral (4). 

From the study of the periodic orbits starting perpen- 
dicularly from the x axis in the above four resonance cases 
we can say the following. In the 1:1 resonance case there 
are analytical results for the position and period of the 
periodic orbits. On the other hand, in the other three res- 
onance cases the semi-numerically derived formulas give 
interesting results. The outcomes from these formulas are 
in very good agreement with the results from numerical 
integration in the cases of the 2:3 and 4:3 resonant cases, 
while in the case of the 2:1 resonance the observed devia- 
tions are smaller than 10%. 



3. Periodic orbits going through the origin 

In this Section we shall study the periodic orbits going 
through the origin. We start again from the he 1:1 inter- 
esting resonant case, when u>i = 0J2 = u). For the motion 
along the lines 



(21) 



the equations of motion become identical, while the fre- 
quencies of the oscillations are always equal. Thus the 
straight lines (21) are exact 1:1 resonant periodic orbits 
going though the origin. The value of the corresponding 
period can be found if we use the formula 



W s 



w s = 



The frequency w s is given by the equation 
, 1/2 

.2 



ui 



3/2 



(22) 



(23) 



where r a is the radius of the circle V(r) = h. Therefore all 
iso-energetic orbits (21) have the same period. In Table 
4 we make a comparison between the period T s and the 
period T„ found using numerical integration, for the 1:1 
straight line periodic orbits, for six values of the energy h. 
As one can see there is a good agreement between T s and 

In the case of the 2:1 resonance we did not find periodic 
orbits going through the origin. In fact, we searched up to 
values of the energy five times as large as those of Table 2, 
without finding any 2:1 periodic orbits going through the 
origin. 

Thus we come to find the periodic orbits going through 
the origin in the 2:3 resonance. The numerical calculations 
support the idea that the orbit starts from the origin with 



Pxs = \ — [h+ — 
tdi V a 



(24) 



while the value of p ys is found from the energy integral 
(4). At the origin the system behaves as separable. The 
energy E xs along the x axis is 



V 



(25) 



4 



Tabic 5: Positions and periods of the 2:3 resonant periodic or- 
bits going through the origin. The values of the parameters are: 
lui = 0.4, u>2 = 0.6. Subscript s indicates values derived by semi- 
numerical methods, while subscript n shows results obtained by nu- 
merical integration. 



h 


Pxs 


Pxn 


Vs 


Vn 


T s 


T n 


0.065 


0.3217 


0.3319 


0.2713 


0.2652 


31.2222 


31.1287 


0.070 


0.3331 


0.3421 


0.2838 


0.2803 


31.2360 


31.1513 


0.075 


0.3442 


0.3520 


0.2958 


0.2947 


31.2483 


31.1711 


0.080 


0.3549 


0.3618 


0.3073 


0.3079 


31.2592 


31.1885 


0.085 


0.3653 


0.3672 


0.3184 


0.3331 


31.2690 


31.2016 


0.090 


0.3755 


0.3768 


0.3291 


0.3444 


31.2778 


31.2157 



Thus 



E, 



ys 



h-E x 



(26) 



We also observe that these orbits intersect the y axis 
perpendicularly. Let y s be the point of intersection. Nu- 
merical experiments show that the point of intersection 
can be found with sufficient accuracy using the equation 



Vs 



2h 



P 2 



1/2 



(27) 



The period of the 2:3 periodic orbits going through the 
origin is well approximated if we use the formula 



T„ = 



An 



10(y2 + Q 2)3/2 



1/2- 



(28) 



Table 5 gives the values of the starting positions p xs 
and p xn of the 2:3 periodic orbits going through the ori- 
gin. In the same Table, we give the values of y s and y n 
and the corresponding periods T s and T n for six values of 
the energy h. The agreement is good. The values of the 
parameters and the energies are as in Table 3. Figure 4 
shows a 2:3 periodic orbit going through the origin for the 
above values of parameters when h = 0.09. Initial condi- 
tions are: x = yo = 0,p xn = 0.3768, while the value of 
p y o is found from the energy integral (4). 

Our last but not least case, is to give semi-numerical 
results for the position of the 4:3 periodic orbits going 
through the origin. It is found that the starting position 
of these orbits is given with good agreement if we use the 
following value for the p x 



Px 



= , l—h, 



(29) 



while the value of p ys is found from the energy integral (4). 
Working as in the case of the 2:3 resonance we find that 
the y s can be obtained by the same formula (27), while a 
good approximation for the period is obtained if we use 
the equation 



8tt 



2(3/2 +a 2) 



3/2 



1/2- 



(30) 
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Figure 4: A 2:3 periodic orbit going through the origin, when h = 
0.09. The values of parameters arc as in Table 2. Initial conditions 
arc: xq = yo = 0,p x o = p X n = 0.3768, while the value of p y o is 
found from the energy integral. 



Table 6: Similar to Table 5 for the 4 
values of the parameters arc: uj\ = 



3 resonant periodic orbits. The 
8,w 2 = 0.6. 



h 


Pxs 


Pxn 


Vs 


Vn 


T s 


T 


0.35 


0.5123 


0.4975 


1.1024 


1.1232 


31.4074 


31.3698 


0.40 


0.5477 


0.5373 


1.1785 


1.1934 


31.4089 


31.3775 


0.45 


0.5809 


0.5725 


1.2500 


1.2624 


31.4100 


31.3830 


0.50 


0.6123 


0.6021 


1.3176 


1.3317 


31.4108 


31.3873 


0.55 


0.6422 


0.6326 


1.3819 


1.3957 


31.4115 


31.3907 


0.60 


0.6708 


0.6624 


1.4433 


1.4553 


31.4120 


31.3935 



Table 6 is similar to Table 5 for the 4:3 periodic orbits 
going through the origin. The values of the parameters 
and the energies are as in Table 3. Figure 5 shows a 4:3 
periodic orbit going through the origin for the above values 
of parameters when h = 0.50. Initial conditions are: xq = 
yo = 0,p xn = 0.6021, while the value of p y o is found from 
the energy integral (4). 

In this Section we have studied the orbits going through 
the origin. Here we presented the resonant cases 1:1, 2:3 
and 4:3. In the resonant case 2:1 we did not find reso- 
nant periodic orbits going through the origin. Our study 
shows that the semi-numerical formulas can give good a 
approximation for the position and periods of the resonant 
periodic orbits. 

4. Discussion 

In this paper we have studied the periodic orbits in 
a dynamical system composed of a two dimensional har- 
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Figure 5: A 4:3 periodic orbit going through the origin, when h = 
0.50. The values of parameters are as in Table 3. Initial conditions 
arc: xq = yo = 0,p x o = p X n = 0.6021, while the value of p y o is 
found from the energy integral. 

monic oscillator and a Plummer potential in the resonant 
cases u>i : UJ2 — n/m = 1:1, 2:1, 2:3 and 4:3. We chose 
the above resonances because all of them are important 
resonances. On the other hand, it is not possible to study 
all resonances in a given potential. 

In this potential we presented semi-numerical formulas 
for the position and the period of the resonant periodic 
orbits intersecting the x axis perpendicularly and going 
through the origin. In semi-numerical methods we com- 
bine empirical or theoretical results or both, with results 
found numerically, in order to find relations or formulas 
giving outcomes that are close to those obtained by pure 
numerical methods, that is by numerical integration. At 
this point, it would be useful to refer to how the numerical 
results indicate or suggest the formulas giving the position 
and period of the periodic orbits. In other words, to ex- 
plain how the formulas are obtained. This is done on a 
complete empirical basis, using simple expressions which, 
of course, contain the parameters entering the Hamilto- 
nian (4) and give outcomes that are close to those given 
by the numerical integration. If the agreement is not sat- 
isfactory, we look for a new formula and so on, until we 
have achieved a good result. In some cases, we combine 
empirical and theoretical results in order to obtain a better 
accuracy (see equations (11)-(14)). 

We have presented semi-numerical results in previous 
papers (see Paper I; Caranicolas, 1994a, b) not only for 
polynomial potentials representing local galactic motion 
(i.e. near an equilibrium point) but also for global galactic 
motion (Caranicolas and Innancn, 1991). Semi-numerical 



methods were also used in Celestial Mechanics (see Hcn- 
rard and Caranicolas, 1990). Therefore these methods are 
a sharp tool in the study of dynamical systems. 

It is well known that for the study of dynamical sys- 
tems, scientists have used numerical or analytical methods 
(Deprit, 1991; Deprit and Elipe, 1991; Elipe, 2000; Elipe 
and Deprit, 1999), or the integrals of motion (Ferrer et 
al., 1997; Ferrer et al., 1998; Ferrer et al., 1998) or maps 
(Caranicolas, 1990). But analytical methods are not easy 
to apply in non-polynomial potentials (i.e. global galactic 
potentials such as mass models). Moreover maps or for- 
mal integrals of motion are not easy to find for potentials 
of the form (1). Therefore we decided to apply the semi- 
numerical method in the potential (1). The good agree- 
ment between the results of the semi-numerical method 
and the outcomes given by the numerical integration of 
the equations of motion justifies our choice. 
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